PROGRAM test9j

USE rrf_module
USE wigner
USE input
IMPLICIT NONE


INTEGER :: a, b, c, d, e, f, g, h, i, m, n, p, x, minj=0, maxj=12,      &
    rule, count(0:1)=0
TYPE(rrf), POINTER :: k=>null(), k1=>null(), k2=>null(), k3=>null(),    &
    k4=>null(), sum=>null(), z=>null(), unity=>null()
LOGICAL :: eof, verbose=.false., printall=.false., allok, single=.false.
CHARACTER(LEN=16) :: key

write (0,"(a)") "Please wait ..."

call init_rrf(.true.,100000)
call unit(unity)
! print "(i0)", unity%factor

data: do
  call read_line(eof)
  if (eof) exit
  call reada(key)
  select case(upcase(key))
  case default
    call die("Keyword "//trim(key)//" not recognized",.true.)
  case("","!","NOTE")
    cycle
  case("VERBOSE")
    verbose=.true.
  case("QUIET")
    verbose=.false.
  case ("PRINT")
    call readu(key)
    select case(key)
    case("ALL")
      printall=.true.
    case("FAILED","ERRORS")
      printall=.false.
    end select
  case("MIN")
    call readi(minj)
    minj=2*minj
  case("MAX")
    call readi(maxj)
    maxj=2*maxj

  case("SUM-RULE","SUMRULE","SUM", "TEST")
    rule=1
    single=.false.
    do while (item<nitems)
      call readu(key)
      select case(key)
      case("RULE")
        call readi(rule)
      case("SINGLE")
        single=.true.
      case default
        call reread(-1)
        call readi(rule)
      end select
    end do
    select case(rule)

    case(1)
      !  Sum rule 1 (orthogonality)

      !  sum_{cf} (2c+1)(2f+1)(2g+1)(2h+1) {a b c}{a b c}
      !                                    {d e f}{d e f}
      !                                    {g h i}{m n i}
      !   = delta_{gm} delta{hn}

      x=2
      allok=.true.
      print "(/a)", "Sum rule 1"
      if (printall .or. single)                                       &
          print "(/a)", " a    b    d    e    g    h    i    m    n"
      if (single) then
        do
          call read_line(eof)
          call reada(key)
          if (upcase(key) == "END") then
            exit
          else
            call reread(1)
          end if
          if (eof) exit data
          if (nitems==0) cycle
          call readargs(a,b,d,e,g,h,i,m,n)
          call sumrule1(a,b, d,e, g,h, i, m,n, .true.)
        end do
      else
        count=0
        do a=minj,maxj
          do b=minj,a
            do d=minj,a
              do e=minj,a
                do g=max(minj,abs(a-d)),min(a+d,maxj),2
                  do h=max(minj,abs(b-e)),min(b+e,maxj),2
                    do i=max(minj,abs(g-h)),min(g+h,maxj),2
                      do m=max(minj,abs(a-d)),min(a+d,maxj),2
                        do n=max(minj,abs(b-e)),min(b+e,maxj),2
                          call sumrule1(a,b, d,e, g,h, i, m,n, printall)
                        end do
                      end do
                    end do
                  end do
                end do
              end do
            end do
          end do
        end do
      end if
      print "(a,i0,a,i0)", "Successful ", count(0), ",  failed ", count(1)
      if (allok) print "(a)", "Test completed -- no errors"
      call clearall

    case(2)
      !  Sum rule 2 (contraction of 6j symbols)

      !  {a b c} = sum_{k} (-)^{2k}(2k+1) {a b c}{d e f}{g h i}
      !  {d e f}                          {f i k}{b k h}{k a d}
      !  {g h i}

      x=2
      allok=.true.
      print "(/a)", "Sum rule 2"
      if (printall .or. single)                                       &
          print "(/a)", " a    b    c    d    e    f    g    h    i"
      if (single) then
        do
          call read_line(eof)
          call reada(key)
          if (upcase(key) == "END") then
            exit
          else
            call reread(1)
          end if
          if (eof) exit data
          if (nitems==0) cycle
          call readargs(a,b,c,d,e,f,g,h,i)
          call sumrule2(a,b,c,d,e,f,g,h,i, .true.)
        end do
      else
        count=0
        do a=minj,maxj
          do b=minj,a
            do d=minj,a
              do e=minj,a
                do c=max(minj,abs(a-d)),min(a+d,maxj),2
                  do f=max(minj,abs(d-e)),min(d+e,maxj),2
                    do g=max(minj,abs(a-d)),min(a+d,maxj),2
                      do h=max(minj,abs(b-e)),min(b+e,maxj),2
                        do i=max(minj,abs(g-h)),min(g+h,maxj),2
                          call sumrule2(a,b,c,d,e,f,g,h,i,printall)
                        end do
                      end do
                    end do
                  end do
                end do
              end do
            end do
          end do
        end do
      end if
      print "(a,i0,a,i0)", "Successful ", count(0), ",  failed ", count(1)
      if (allok) print "(a)", "Test completed -- no errors"
      call clearall

    end select
  end select
end do data


CONTAINS

SUBROUTINE sumrule1(a,b, d,e, g,h, i, m,n, print)

INTEGER, INTENT(IN) :: a,b,d,e,g,h,i,m,n
LOGICAL, INTENT(IN) :: print

INTEGER :: c, f, x
LOGICAL ok
CHARACTER(LEN=45) :: args
CHARACTER(LEN=8) :: verdict
CHARACTER(LEN=100) :: string

call unit(unity)
! call list(unity)
x=2
call setzero(sum)
do c=abs(a-b),a+b,2
  do f=max(abs(d-e),abs(c-i)),min(d+e,c+i),2
    call ninej(a, b, c, d, e, f, g, h, i, x, k1)
    call ninej(a, b, c, d, e, f, m, n, i, x, k2)
    call int_to_rrf(int((c+1)*(f+1)*(g+1)*(h+1),long),k)
    call mult(k,k1)
    call mult(k,k2)
    call add(sum,k)
  end do
end do
! call list(sum)
! call copy(sum,k)
if (g==m .and. h==n) then
  ok=equal(sum,unity)
!  call subtr(k,unity)
else
  ok=iszero(sum)
end if
! ok=iszero(k)
if (ok) then
  verdict="  O.K."
  count(0)=count(0)+1
else
  verdict=" Wrong!"
  if (allok .and. .not. printall)            &
      print "(/a)", " a    b    d    e    g    h    i    m    n"
  allok=.false.
  count(1)=count(1)+1
end if
if (print .or. .not. ok) then
  p=1
  args=""
  call put(a, args, p)
  call put(b, args, p)
  call put(d, args, p)
  call put(e, args, p)
  call put(g, args, p)
  call put(h, args, p)
  call put(i, args, p)
  call put(m, args, p)
  call put(n, args, p)
  p=1
  string=""
  call char4i(sum, string, p, 40)
  print "(4a)", args, verdict, "  Sum = ", trim(string)
end if

END SUBROUTINE sumrule1

!-----------------------------------------------------------------------

SUBROUTINE sumrule2(a,b,c,d,e,f,g,h,i, print)

INTEGER, INTENT(IN) :: a,b,c,d,e,f,g,h,i
LOGICAL, INTENT(IN) :: print

INTEGER :: j, x
LOGICAL ok
CHARACTER(LEN=45) :: args
CHARACTER(LEN=8) :: verdict
CHARACTER(LEN=100) :: string

call unit(unity)
! call list(unity)
x=2
call setzero(sum)
do j=max(abs(b-f),abs(d-h),abs(a-i)),min(b+f,d+h,a+i),2
  call sixj(a, b, c, f, i, j, x, k1)
  call sixj(d, e, f, b, j, h, x, k2)
  call sixj(g, h, i, j, a, d, x, k3)
  call int_to_rrf(int((j+1),long),k)
  call mult(k,k1)
  call mult(k,k2)
  call mult(k,k3)
  if (mod(j,2) > 0) call chsign(k)
  call add(sum,k)
end do
! call list(sum)
! call copy(sum,k)
call ninej(a,b,c,d,e,f,g,h,i, 2, k4)
ok=equal(sum,k4)
if (ok) then
  verdict="  O.K."
  count(0)=count(0)+1
else
  verdict=" Wrong!"
  if (allok .and. .not. printall)            &
      print "(/a)", " a    b   c   d    e   f    g    h    i"
  allok=.false.
  count(1)=count(1)+1
end if
if (print .or. .not. ok) then
  p=1
  args=""
  call put(a, args, p)
  call put(b, args, p)
  call put(c, args, p)
  call put(d, args, p)
  call put(e, args, p)
  call put(f, args, p)
  call put(g, args, p)
  call put(h, args, p)
  call put(i, args, p)
  p=1
  string=""
  call char4i(sum, string, p, 100)
  string(25:)="  9j = "
  p=32
  call char4i(k4, string, p, 100)
  print "(4a)", args, verdict, "  Sum = ", trim(string) 
end if

END SUBROUTINE sumrule2

!-----------------------------------------------------------------------

SUBROUTINE readargs(a,b,d,e,g,h,i,m,n)

INTEGER, INTENT(OUT) :: a,b,d,e,g,h,i,m,n

call readarg(a)
call readarg(b)
call readarg(d)
call readarg(e)
call readarg(g)
call readarg(h)
call readarg(i)
call readarg(m)
call readarg(n)

END SUBROUTINE readargs

!-----------------------------------------------------------------------

SUBROUTINE readarg(a)

INTEGER, INTENT(OUT) :: a
INTEGER :: p, x
CHARACTER(LEN=10) :: buff

call reada(buff)
p=1
x=2
call number(buff,p,a,x)
if (x==1) a=2*a

END SUBROUTINE readarg

!-----------------------------------------------------------------------

SUBROUTINE put(j, s,p)

CHARACTER(LEN=*) :: s
INTEGER, INTENT(IN) :: j
INTEGER, INTENT(INOUT) :: p

if (mod(j,2) == 0) then
  write (unit=s(p:p+1),fmt="(i2)") j/2
else
  write (unit=s(p:p+3),fmt="(i2,a2)") j,"/2"
endif
p=p+5

END SUBROUTINE put

!-----------------------------------------------------------------------

SUBROUTINE clearall

call clear(k)
call clear(k1)
call clear(k2)
call clear(sum)
call clear(z)
call clear(unity)

END SUBROUTINE clearall

END PROGRAM test9j
